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Abstract. This paper is concerned with the problem of shape optimization of two- 
dimensional flows governed by the time-dependent Navier-Stokes equations. We derive 
the structures of shape gradients with respect to the shape of the variable domain for 
time-dependent cost functionals by using the state derivative with respect to the shape 
of the fluid domain and its associated adjoint state. Finally we apply a gradient type 
algorithm to our problem and numerical examples show that our theory is useful for 
practical purpose and the proposed algorithm is feasible in low Reynolds number flow. 
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1 Introduction 

The problem of finding the optimal design of a system governed by the incompress- 
ible Navier-Stokes equations arises in many design problems in aerospace, automotive, 
hydraulic, ocean, structural, and wind engineering. Example applications include aero- 
dynamic design of automotive vehicles, trains, low speed aircraft, sails, and hydrody- 
namic design of ship hulls, turbomachinery, and offshore structures. In many cases, the 
flow equations do not admit steady-state solutions, and the optimization model must 
incorporate the time-dependent form of the Navier-Stokes equations. 

Optimal shape design has received considerable attention already. Early works 
concerning on existence of solutions and differentiability of the quantity (such as, state, 
cost functional, etc.) with respect to shape deformation occupied most of the 1980s (see 
[2, 3, 15, 16, 21]), the stabilization of structures using boundary variation technique 
has been fully addressed in [3, 16, 21]. However, a few studies have considered the 
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shape optimization of time-dependent flows (see [5, 11, 14, 19]). Our concern in this 
article is on shape sensitivity analysis of time-dependent Navier-Stokes flow with small 
regularity data, and on deriving an efficient numerical approach for the solution of 
two-dimensional realizations of such problems. 

In [7], we use the state derivative approach to solve a shape optimization problem 
governed by a Robin problem, and in [8, 9], we derive the expression of shape gradients 
for Stokes and Navier-Stokes optimization problem by this approach, respectively. In 
this paper, we use this approach and weak implicit function theorem to derive the 
structures of shape gradients with respect to the shape of the variable domain for 
some given cost functionals in shape optimization problems for time-dependent Navier- 
Stokes flow with small regularity data. 

This paper is organized as follows. In section 2, we briefly recall the velocity method 
which is used for the characterization of the deformation of the shape of the domain 
and give the definitions of Eulerian derivative and shape derivative. We also give the 
description of the shape optimization problem for the time-dependent Navier-Stokes 
flow. 

In section 3, we employ the weak implicit function theorem to prove the existence of 
the weak Piola material derivative, and then give the description of the shape deriva- 
tive. After that, we express the shape gradients of some typical cost functionals by 
introducing the corresponding linear adjoint state systems. 

Finally in section 4, we propose a gradient type algorithm with some numerical 
examples to prove that our theory could be very useful for the practical purpose and 
the proposed algorithm is efficient in low Reynolds number flow. 

2 Preliminaries and statement of the problem 

2.1 Elements of the velocity method and notations 

Domains don't belong to a vector space and this requires the development of shape 
calculus to make sense of a "derivative" or a "gradient". To realize it, there are about 
three types of techniques: J.Hadamard [10] 's normal variation method, the perturbation 
of the identity method by J.Simon [17] and the velocity method (see J.Cea[2] and J.- 
P.Zolesio[3, 20]). We will use the velocity method which contains the others. In that 
purpose, we choose an open set D in with the boundary dD piecewise C^, and a 
velocity space V eE'' ■.= {V e C{[0,e];V''{D,R^)) : V-uqd = on dD}, where e is a 
small positive real number and D'^(I),M^) denotes the space of all /c— times continuous 
differentiable functions with compact support contained in . The velocity field 

V{s){x) = V{s,x), x£D, s>0 

belongs to V''{D,R^) for each s. It can generate transformations 

Ts{V)X = x{s,X), s>0, XeD 
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through the following dynamical system 

(^{s,X) = V{s,x{s)) 
\xiO,X)=X 

with the initial value X given. We denote the "transformed domain" Ts(V){^l) by 
fls(V) at s > 0, and also set dQg ■= Ts{di}). 

There exists an interval I = [0,6), < 6 < e, and a one-to-one map Tg from D onto 
D such that 

(i) To = I; 

(ii) {s,x) ^ Ts{x) belongs to C^{I;C^{D;D)) with Ts{dD) = dD; 

(iii) {s,x) ^ r-i(x) belongs to C{r,C^{D; D)). 

Such transformation are well studied in [3]. 

Furthermore, for sufficiently small s > 0, the Jacobian Jg is strictly positive: 

Jsix) := det|DT,(2;)| = detDT,(x) > 0, (2.2) 

where DT<j(x) denotes the Jacobian matrix of the transformation Tg evaluated at a point 
x (z D associated with the velocity field V. We will also use the following notation: 
DT~^(x) is the inverse of the matrix DTs{x) , *DT~^(x) is the transpose of the matrix 
Dr~^(x). These quantities also satisfy the following lemmas. 

Lemma 2.1 ([21]) For any V ^ , DT^ and Js are invertible. Moreover, DTg, 
DT-i are m Ci([0, e]; C'=-i(£'; ]R^><^)), andJs, are in C\[0,e];C''-^D;R)) 

Lemma 2.2 ([21]) cp is assumed to be a vector function in C^{D)^ . 

(1) D(T7i)oT, = DT-i; 

(2) D((^ o r-i) = ■ DT-i) o T-i; 

(3) (D</.)or, = D((^oT,)-DT-i. 

Now let J(0) be a real valued functional associated with any regular domain 17, we 
say that this functional has a Eulerian derivative at in the direction V if the limit 

lim ■^<»'> - ■^<»' := dJ(0;V) 

s\0 S 

exists. 

Furthermore, if the map 



V ^ dJ{n-V) : E'' 
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is linear and continuous, we say that J is shape differentiable at Vl. In the distribu- 
tional sense we have 



dJ(17;l/) = (VJ,F)^.(5,; 



(2.3) 



When J has a Eulerian derivative, we say that VJ is the shape gradient of J at VL. 

Before closing this subsection, we introduce the following functional spaces which 
will be used throughout this paper: 

if(div,17) :={MeL2(17)^ : div'a = OinJl, -a • n = on 917}, 
ifo^(div,17) := {it G H^(yt)^ : divit = in u\dn = 0}. 

Given T > 0, we introduce the notation LP{0,T; X) which denotes the space of 
integrable functions / from [0, T] into the Banach space X with the norm 

\lp{o,T;X) = ( [ 11/11x^*1 ' l<P<+oo. 



We also denote by L°°(0, T; X) the space of essentially bounded functions / from [0, T] 
into X, and is equipped with the Banach norm 

ess sup 

te[o,T] 

2.2 Statement of the shape optimization problem 

In two dimensions, we consider a typical problem in which a solid body S with the 
boundary dS is located in an external flow. Since the flow is in an unbounded domain, 
we reduce the problem to a bounded domain D by introducing an artificial boundary 
dD on which we set the speed flow y = y^. := D\S is the effective domain with 
its boundary d^l = dS U dD. The state equations of the flow can be written by the 
Navier-Stokes equations in the non-dimensional form. 



(2.4) 



dty - aAy + Dy ■ y + Vp = 
div y = 

y = yoo 

y = 

y(o) = vo 

f^pdx = 0, 
JoDVoo-nds = 0, 

where the last relation is needed in view of the incompressibility constraint divy = 0, a 
stands for the inverse of the Reynolds number whenever the variables are appropriately 
nondimensionalized, y,p, and / are the velocity, pressure, and the given body force per 
unit mass, respectively. 

Our goal is to optimize the shape of the boundary dS which minimizes a given cost 
functional J depending on the fluid state. The cost functional may represent a given 



in Q := 17 X (0,r), 
in Q, 

on dD X (0,r), 
ondS X (0,r), 
in 17, 

on (o,r), 
on (o,r). 
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objective related to specific characteristic features of the fluid flow (e.g., the deviation 
with respect to a given target velocity, the drag, the vorticity, ...). 

Hence, we are interested in solving the following minimization problem 

mmJi{n) = \[ [\y-yfdxdt, (2.5) 



or 



mmJ2{n) = - [ /icurlypdxdt, (2.6) 

where y is satisfied by the full Navier-Stokes system (2.4) and y^ is the target velocity 
given by the engineers. We also notice that the boundary dD is fixed in our optimization 
problems and an example of the admissible set O is: 

O := i C : dD is fixed, dx = constant 



In order to deal with the nonhomogeneous Dirichlet boundary condition on dD, let 
the vectorial function h be the solution of 

div h = in 17 

h = y^ on dD (2.7) 
h = ondS, 
then we can choose an extension h with h, = in the body S. 

Now we may look for a solution of the nonhomogeneous Navier-Stokes equations in 
the form 

y = h + y, (2.8) 

with y vanishing on the boundary of the domain il. Substituting (2.8) in the system 
(2.4), we find the following equations for y: 

dty - aAy + Dy ■ y + I)y ■ h + Bh ■ y + Vp = F in Q, 
divt/ = inQ, 
y = ondnx (0,r), 

^y{o) = yo inn 

where F := f + a^h — Dh ■ h and := y^ — h. 

For the existence and uniqueness of the solution of the full Navier-Stokes system 
(2.9), we have the following results (see [18]). 

Theorem 2.1 The domain $7 is supposed to be piecewise . We assume that 

f,dtf eL''{0,T;H{diY,D)), (2.10) 

yoGHHD)^nH^{div,D), (2.11) 

y^ G Hy^dD)^, (2.12) 

and the solution of (2.4) is unique and satisfies 

y, dty G L\0, T; ( div , 17)) n L°°(0, T; H{ div , J7)). 

Moreover, if Q is of class and f € L°°{0,T; H^div , D)), then the function y € 
L'^{0,T;H^{Q)^). 
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3 State derivative approach 

In this section, we shall prove the main theorem using an approach based on the differ- 
entiability of the solution of the Navier-Stokes system (2.9) with respect to the variable 
domain. To begin with, we use the Piola transformation to bypass the divergence free 
condition and then derive a weak material derivative by the weak implicit function the- 
orem. Then we will derive the structure of the shape gradients of the cost functionals 
by introducing the associated adjoint state equations. 



3.1 Piola material derivative 

From now on, we assume that n is of class and (2.10)-(2.12) hold. Then we say 
that the function y G L'^{0, T; i^Q ( div , Q)) is called a weak solution of problem (2.9) if 
it satisfies 

{e{y),w) = 0, weL\0,T;H^{div,n)), (3.1) 
with e{y) := {ei{y),e2iy)), := (0,0), and 

{ei{y),w) := / / {dty-w + al)y :Dw+Dy-y-w+Dy-h-w + Dh-y-w — F-w)dxdt, 
Jo Jn 

(3.2) 

{e2{y),w):= [ {y{0)-yo)MO)dx. (3.3) 
Jn 

It must be considered that the divergence free condition is variant with respect to the 
use of the transformation Ts during the derivation of the shape gradient for the cost 
functional. Therefore, we need to introduce the well known Piola transformation which 
preserves the divergence free condition. 

Lemma 3.1 ([1]) The Piola transform 

: H{div ,n)^H{div 

is an isomorphism. 

Now by the transformation Tg, we consider the solution y^ defined on x (0, T) 
of the perturbed weak formulation 

/ / {dtys ■ Ws + aDj/g : Bws + Dy^ • • lu^ + By^ h ws 

Jo JQs 

+ 'Dh-y,-Ws- F ■Ws)dxdt = 0, (3.4) 

/ {yM-yo)-'Ws{0)dx = 0, (3.5) 
JQs 
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for a\l Ws G ^^(0, T; i^g ( div , O^)), and introduce = ^'^^(j/J, w"* = ^'7^(10^) defined 
on Q. Then we replace y^, Wg by ^'^(j/*), ^'^(tu*) in the weak system (3.4)(3.5) 

^ [ m^sir) ■ ^.K) + «D(^,(r)) : Bi^siw)) 

+D(^,(r)) • ^s(r) • ^.(«^') + D(^.(r)) • h ■ ^s{w') 

+Dh ■ • "ifsiW) - F ■ "^siw')] dxdt = 0, (3.6) 

/ [^,(r(0))-yo]-^.(«^'(0))dx = (3.7) 

Jfls 

for all e L^iO,T; H^{div 

Using a back transport into 0, and employing Lemma 2.2, we obtain the following 
weak formulation 

{e{s,y'),w') = 0, yw' e L'^{0,T-H^{div,n)) (3.8) 

with the notations e := (61,62), where 
rT 



{ei{s,v),w) := [ I dt{B{s)v) ■ {DTsw)dxdt 
Jo Jn 

+ a / B{B{s)v) :[B{B{s)w) ■ A{s)]dxdt+ / / B{B{s)v) ■ v ■ {B{s)w) dx dt 
Jo Jn Jo Jn 

+ [ [ B{B{s)v) -h- {B{s)w)dxdt+ [ [ B{B{s)h) -v {B{s)w)dxdt 
Jo Jn Jo Jn 

-[ [ {FoTs)-(DTs-w)dxdt, (3.9) 
JO Jn 

and 

(62(5, v),w) := / iB{s)v{0) - y^ o T,) ■ (Dr,ty(0)) dx (3.10) 
Jn 

and 

A{s) := JsBT-^*BT-^; B{s)t := J^^DT, • r. 
Now we are interested in the differentiability of the mapping 

s^r = ^;Hys) ■■ [0,e]^L\0,T;H^{dw,n)) 

where e > is sufficiently small and y^ is the solution of the weak formulation 

{e{s,v),w) = 0, yw e L^iO,T-H^{div ,n)). (3.11) 

In order to prove the differentiability of with respect to s in a neighborhood of s = 0, 
there maybe two approaches: 

(i) analysis of the differential quotient: lim(j/* — y)/s] 

s^O 



(ii) derivation of the local differentiability of the solution y associated to the implicit 
equation (3.1). 

We use the second approach. Since / e L2(0, T; H{ div , D)), we deduce that (/ o - 
f)/s weakly converges to D/ • V in L^(0, T; H^^{D)^) as s goes to zero. Thus we can 
not use the classical implicit function theorem, since it requires strong differentiability 
results in H~^. Hence we introduce the following weak implicit function theorem. 

Theorem 3.1 ([20]) Let X, Y' be two Banach spaces, I an open bounded set in M, 
and consider the map 

is,x) 1-^ e{s,x) : I x X Y' 
If the following hypothesis hold: 

(i) s I— > {e{s,x),y) is continuously differentiable for any y £ Y and{s,x) i— > {dse{s , x) , y) 
is continuous; 

(ii) there exists u £ X such that u G C^'^(r,X) and e(s,n(s)) = 0, Vs G I; 
(Hi) X ^ e{s,x) is differentiable and {s,x) i— > dxe{s,x) is continuous; 

(iv) there exists sq € / such that dxe{s, x)\(^g^^x{so)) 'i-^ o,n isomorphism from X to Y' , 
the mapping 

s I— > u{s) : / I— > X 

is differentiable at s = sq for the weak topology in X and its weak derivative u{s) is the 
solution of 

{dxe{so,uiso)) ■ u{so),y) + {dse{so,u{so)),y) = 0, Vy G Y. 
We may now state the main theorem of this section. 

Theorem 3.2 We assume that the domain is piecewise and (2.10) -(2.12) hold, 
y G L^(0, T; //q ( div , ri)) is the solution of the weak formulation (3.1). Then the weak 
Piola material derivative y^ := ds{y^)\s=o exists and is characterized by the follow- 
ing weak formulation: 

{dveiO,v)\^^y ■ y^,w) + {dse{Q,y),w) = 0, e L\Q,T- Hl{diY ,^)), (3.12) 
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I.e., 



I / [Sf^^-w + aDy^ : Dw + Dy^-j/-iu + Djif-y^-i(;+Dy^-/i-itJ + D/i-^^-iu] dxdt 
Jo Jn 

= - I I [dt{(DV - div V)y) -w + dty -BV ■w]dxdt 
Jo Jn 

-a / B{(DV - divV)y) :Dwdxdt 
Jo Jn 

-a I I By : [B{{BV - div V)w) +Bw ■ {divV - BV - *BV)] dx dt 
Jo Jn 

- [ [ [D((DV- divV){y + h)) ■ y ■ w + B{y + h) ■ y ■ {{BV - divV)w)]dxdt 
Jo Jn 

rp 

- I I [D((DF - div V)y) hw-Byh - {{BV - div V)w)] dx dt 
Jo Jn 

+ I I {*BV ■ (/ + aAh -Bh h) + B{f + aAh - Bh ■ h) ■ V) ■ w dx dt 
Jo Jn 

+ I if + -Bhh)- {BV ■ w) dx dt, (3.13) 
Jo Jn 

and 

[ ii^{0)-w{0)dx = - [ [{BV + *BV - divVl)-y{0)-{ByQV + *BVyQ)]-w{0)dx. 
Jn Jn 

(3.14) 

Proof. In order to apply Theorem 3.1, we need to verify the four hypothesis of 
Theorem 3.1 for the mapping 

{s,v) ^ e{s,v) : [0,e] x L'^{0,T; H^{dw ^ L'^{0,T; H^{dw ,Qy). 

To begin with, since Q is of piecewise C^, the mapping Ts € C^([0, e]; C^{D, D)). Then 
by Lemma 2.1, the mapping 

s (ei(s,-u),w) : [0,e] M (^ = 1,2) 

is for any v,w ^ L2(0, T; i?|J( div , $7)). On the other hand, since / G L?{'d,T; 
H{dw , D)), the mapping s i— > / o is only weakly differ entiable in H~^, thus the 
mapping s ei{s,v) is weakly differ entiable, and then s i— > e{s,v) is weakly differen- 
tiable. 

Since we have the following identities by simple calculation, 

-^DT, = (Dy(s)or,)Dr,; (3.15) 
ds 

-^Js = {diyV{s))oT,Js; (3.16) 
ds 

^(/or,) = (D/. v(s))or„ (3.17) 
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the weak derivative of ei{s, v) (i = 1, 2) can be expressed as 

{dsei{s,v),w) = r [ [dt{B'{s)v)-{BTsw) + dt{B{s)v)-iBV{s)oT,)-BTs-w]dxdt 
Jo Jn 

+ a [ [ B{B'{s)v) : [DiB{s)w) ■ A{s)]dxdt 
Jo Jn 

+ a [ [ D{B{s)v) : [B{B'{s)w) ■ A{s) + I){B{s)w) ■ A'{s)] dxdt 
Jo Jn 

+ I I [D{B'{s)v) ■ V ■ {B{s)w) + B{B{s)v) ■ v ■ {B'{s)w)] dxdt 
Jo Jn 

+ [ [ [B{B'{s)v) ■ h ■ iB{s)w) + B{B{s)v) ■ h ■ iB'{s)w)] dx dt 
Jo Jn 

+ [ [ [B{B'{s)h) ■ V ■ {B{s)w) + B{B{s)h) ■ v ■ {B'{s)w)] dx dt 
Jo Jn 

- [ [ [*DV{s) -F + BF- V{s)] o Ts ■ (DT, w) dx dt 
Jo Jn 

- r f {FoT,)-[{BV{s)oTs)-BTs-w]dxdt, (3.18) 
Jo Jn 

and 

{dse2{s,v),w) = [ {[B'{s)v{0) - (Dj/o • V{s))oTs] ■ (DT,to(0)) 
Jn 

+ iB{s)v{0) - yo o T,) ■ [{BVis) o T,) • DT, • w{0)]}dx, (3.19) 

where 

B'{s)t := ^ [B{s)t] = [BV{s) o - ( div V{s) o Ts)1]B{s)t; 
A'{s) := f^Ais) = [divF(s) oT, - BT-^BVi-s) oT,] A{s) - *[BT-^BV{s) oT, A{s)]. 

Obviously, the mapping {s,v) i-^ dse{s,v) is continuous, and when we take s = 0, we 
have 

B'{0)t = {BV - divFI)-T; 
^'(0) = div Vl-BV - *BV, 
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and then 

{dsei{0,v),w) = [ [ [dt{(DV - dwV)v) -w + dtv -BV -wjdxdt 
Jo Jn 

+ Q / / BiCDV - dwV)v) -.Bwdxdt 
Jo Jn 

+ a [ [ Bv: [B{(BV - div V)w) + Bw ■ {dWV - BV - *BV)] dx dt 
Jo Jn 

+ [ [ [B{{BV - divV){v + h)) ■ V ■ w + B{v + h) ■ V ■ {{BV - divV)w)]dxdt 
Jo Jn 

+ [ [ [B{{BV - div V)v) h w + Bv h - {{BV - div V)w)] dx dt 
Jo Jn 

-[ [ F ■ {BV ■w)dxdt- [ [ {*BV ■ F + BF -V) -wdxdt. (3.20) 
Jo Jn Jo Jn 

(5,62(0, v),w) = [ [{BV + *BV - div V) ■ v{0) ■ w{0) - {By^V + *BVyo) ■ w{0)] dx. 
Jn 

(3.21) 

To verify (ii), we follow the same steps described in R.Dziri[4] to find that the mapping 
s i/g oTg is Lipschitz continuous which is the direct consequence of the uniqueness 
of the solution of the Navier-Stokes system, i.e., Theorem 2.1. 
It is easy to check that the mappings 

V ^ ei{s,v): L^{0,T;H^{div,n)) L'^{0,T; H^{dw ,ny) 
v^e2{s,v): H^{div ,n) ^ H^{div ,n) 

are differentiable, and the derivatives of ei{s,v) with respect to v in the direction 6v 
are 

rT 



{d^ei{s,v) -Sv^w) = [ [ dt{B{s)6v) ■ {BTsw)dxdt 
Jo Jn 

+ a I I B{B{s)5v) : [B{B{s)w) ■A{s)]dxdt 
Jo Jn 

+ [ [ [B{B{s)5v) ■ V ■ {B{s)w) + B{B{s)v) ■ 6v ■ {B{s)w)] dx dt 
Jo Jn 

+ [ [ [B{B{s)6v)-h-{B{s)w)+B{B{s)h)-5v-{B{s)w)]dxdt. (3.22) 
Jo Jn 



and 



(9^62(5, v) ■ 5v, w)= {B{s)5v{<d)) ■ {BTs ■ w{<d)) dx. (3.23) 
Jn 

The continuity of (s,u) ^ di,ei{s,v) is easy to check. Moreover, 

(9^ei(0, 1?) • (5i7, lu) = / / [dt{5v) ■ w + aB{8v) -.Bw + B{6v) ■ V ■ w 
Jo Jn 

Bv ■ 5v ■ w + B{6v) ■ h ■ w + Bh ■ 6v ■ w]dxdt, (3.24) 
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{d^e2{0, v) ■ 6v,w) = / Sv{0) ■w{0)dx. (3.25) 
Jn 

Furthermore, 6v 8^6(0, v) ■ 5v is an isomorphism which fohows from the uniqueness 
and existence of the Navier-Stokes system, i.e.. Theorem 2.1. Indeed, we assume that 
yi,y2 are two solutions of the Navier-Stokes system (2.9), and y^ (i = 1,2) satisfies 
the weak formulation (3.1). It is obvious that y = yi — y2 satisfies 

/ / [dty w+aT)y -.lyw+By-h-w+Bh y w+By-yi-w+By (3.26) 
JO .m 

and 

/ y(0) ■w{<d)Ax = 0. (3.27) 

Jo. 

Now let = y, we can follow the proof of the unique solvability of the unsteady 
Navier-Stokes equations (see Temam [18]) and obtain 



|y(t)P<o, vte[o,r]. 

Thus yi = 1/2- Similar a priori estimates hold for 5v and the uniqueness of the solution 
of the system (3.24) (3.25) is obtained. 

Finally, all the hypothesis are satisfied by (3.8), we can apply Theorem 3.1 to (3.8) 
and then use (3.20), (3.21), (3.24) and (3.25) to obtain (3.13) and (3.14). □ 

3.2 Shape derivative 

In this subsection, we will characterize the shape derivative y', i.e., the derivative of 
the state y with respect to the shape of the variable domain. 

Theorem 3.3 Under the assumption of Theorem 2.1 and moreover assume that $7 is 
of class C^, y G L°°(0, T; if2(f7)^' n i?,j( div , 17)) solves the weak formulation (3.1) 
and y^ solves the perturbed weak formulation (3.4) (3.5) in Vlg x (0, T), then the shape 
derivative 

y := lim 

exists and is characterized as the solution of 

' dty' - a Ay' + By' ■ y + Dy ■ y' + Dy' ■ h + Dh ■ y' + Vp' = in Q 
div y' = in Q 

y' = -{Dy-n)Vn on dS x {0,T) 

y' = ondDx (0, T) 

y'{0)=0 inn. 

(3.28) 

Proof. Since is of class and V G E^, Jl^ has the same regularity than Q for any 
s E (0,e), then y^ € L°°{0,T; H^{ns)^) satisfies the following weak formulation 

/ / (aDj/g : Bw + By^ ■ygW + By^ ■ h ■ w + Bh ■ y^ ■ w - F ■ w) dx dt = 0, (3.29) 

Jo Jfls 
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/ y^{0) ■w{0)dx = (3.30) 

for any w G L'^{0, T; H^{dw , O^)). Moreover, we have dty, e L'^{0, T; H^{dw , 0)). 
To begin with, we introduce the following Hadamard formula (see [3, 21]) 

/ g{s,x)dx= I ^(s,x)dx+ j g{s,x)V ■ UgdYs, (3.31) 

for a sufficiently smooth functional g : [0, r] x ^ M. 

Now we set a function ip € 'D((5)^~''^ and div(^(x,t) = in for a.e. t G (0,7"). 
Obviously when s is sufficiently small, p{t) belongs to the sobolev space f^Q ( div , ilg) n 
H'^ins)'^ for a.e. t G (0,r). Hence we can use (3.31) to differentiate (3.29), (3.30) with 

[ [ {dtV' + aBy' :B(p + By' ■ y ■ (p + By ■ y' ■ ip + By' ■ h ■ (p + Bh ■ y' ■ if)dxdt 
Jo Jq 

+ / {aBy :Bip + By ■ y ■ If + By ■ h ■ if + Bh ■ y ■ (f - F ■ (fi)Vn ds dt = 0, 
Jo Jan 

[ y'(0) • (^(0) dx+ I yM ■ ^{0)Vn ds = 0. 
Jn Jan 

Since (f has a compact support, the boundary integrals vanish. Using integration by 
parts, we obtain 



[ [ {dty' - a^y' + By' ■ y + By ■ y' + By' ■ h + Bh ■ y') • c^dxdt = 0, 
Jo Jn 



(3.32) 



and 



/ y'(0) • if{0)dx = 0. (3.33) 
Jn 

Then there exists some distribution p' such that 

dty - aAif' + By' ■ y + By ■ y' + By' ■ h + Bh ■ y = -Vp' 

in the distributional sense in Q and y'{0) = in Q since (p{0) is arbitrary. 

Now we recall that for each sufficient small s, ^'^^(yg) belongs to the Sobolev space 
H^{div,n), then we can deduce that its material derivative vanishes on the boundary 
dS. Thus we obtain the shape derivative of y at the boundary dS, 

y' = -ByV, ondSx{0,T). 

Since y\dSx{o,T) = 0, we have D§|g5x(o,T) = D§ • n*n, and then 

y' = -{By-n)Vn ondSx{0,T). 

Since dD is fixed, we obtain y' = on the boundary dD x (0, T). □ 
The shape derivative y' of the solution y of the original Navier-Stokes system (2.4) 
is given by y' = y', then we obtain the following corollary by substituting y' = y' and 
y = y-h into (3.28). 
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Corollary 3.1 The shape derivative y' of the solution y of (2.4-) exists and satisfies 
the following system 

' dty' - aAy' + By' ■ y + By ■ y' + Vp' = in Q; 

diyy' = in Q; 

y' = {-By ■ n)Vn on dS x (0, T) (3.34) 

y' = ondD X (0, T) 

y'{0)=0 inn. 

3.3 Adjoint state system and gradients of the cost functionals 

This subsection is devoted to the computation of the shape gradients for the cost 
functionals Ji{n,) and J2(f^) by the adjoint method. 

For the cost functional Ji(il) = Jq Jq^Iv ~ Vdl'^ dxdt, we have 

Theorem 3.4 Let Q be of class , y^ G L'^{0,T; L^{D)^), and V G the shape 
gradient VJi of the cost functional Ji(0) can be expressed as 



VJi 



-^{y - Vdf + "(Dy • n) ■ {Bv ■ n) 



n. 



(3.35) 



where the adjoint state v satisfies the following linear adjoint system 

—dtV — a/S.v — Bv ■ y + *By ■ v + Vq = y — y^, in Q 

divv = 0, in Q 

V = 0, ondnx (o,r) 

v{T) = 0, in n. 



(3.36) 



Proof. Since Ji(0) is differentiable with respect to y, and the state y is shape 
differentiable with respect to s, i.e., the shape derivative y' exists, we obtain Eulerian 
derivative of Ji^O.) with respect to s, 

dJi(0;V)= r [ iy-ya)-y'dxdt+ T / l\y-yfVndsdt (3.37) 
Jo Jq Jo Jdn ^ 

by Hadamard formula (3.31). 

By Green formula, we have the following identity 

/ / {{dty' - aAy' + By' ■y + Byy' + Vp) ■ w - di\ yV] dx dt 
Jo Jn 

= I I [{—dtw — aAw — Bw ■ y + *By ■ w + Vvr) ■ y' — p' div w] dx dt 
Jo Jn 

+ / {y' ■ w){y ■ n)dsdt+ / {aBw ■ n — Trn) ■ y' ds dt 

Jo Jdn Jo JdQ 

+ [ [ (p'n-aBy'n) ■wdsdt+ [ {y' (T) ■ w{T) - y' (0) ■ w{0)) dx. (3.38) 
Jo Jan Jn 
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Now we define {v,q) to be the solution of (3.36), use (3.34) and set {w,tt) = {v,q) in 
(3.38) to obtain 



/ {y ~ Vd) ■ y' dxdt = — / {aDv ■ n — qn) ■ y ds dt. 

h Jn Jo Jas 



(3.39) 



Since y' = {—Dy ■ Ti)Vn on the boundary dS and divy' = in 0, we obtain the 
Eulerian derivative of Ji(^2) from (3.37), 

rT 



dJi{n;V) 



JdS 



1, 



\y-yd\ + a{Dy ■ n) ■ {Bv ■ n) 



Vn ds dt. 



(3.40) 



Since the mapping V dJi{0,;V) is hnear and continuous, we get the expression 
(3.35) for the shape gradient VJi by (2.3). □ 
For another typical cost functional J2{^) = f /o'" /^l "-^^i have the 

following theorem. 

Theorem 3.5 Let Q, be of class and V G E^, the cost functional J2{^) possesses 
the shape gradient VJ2 which can be expressed as 



VJo = a 



-|curlyp + (Dy ■ n) ■ {Y)v ■ n — curly A n) 



n, 



(3.41) 



where the adjoint state v satisfies the following linear adjoint system 

^ —dtv — aAv — Dv • y + *Dy • v + Vq = —aAy, in Q 

div V = 0, in Q 

v = 0, ondnx{0,T) 

v{T) = 0, in n. 



(3.42) 



Proof. The proof is similar to that of Theorem 3.4. Using Hadamard formula (3.31) 
for the cost functional J2, we obtain the Eulerian derivative 

dJ2i^;V) = a [ [ curly • curly' dxdt+ /" [ -| curlypF„dsdt. (3.43) 
Jo Jn Jo J 80. 2 

Then, we define {v,q) to be the solution of (3.42), use (3.34) and set {w,-k) = {v,q) in 
(3.38) to obtain 

a[ [ Ay -y' dxdt= [ [ a(Dv ■ n) ■ y' ds dt. (3.44) 
Jo Jn Jo Jas 

Applying the following vectorial Green formula 



/ 



{if ■ Atp + curl ip ■ curl i/> + div ip div tp) dx 



/ curl ip An) + (f ■ n div tp) ds 

Jan 



for the vector functions y and y', we obtain 



/ / ( curly • curly' + Ay • y') dx dt = / / (curly A n) • y'dsdt (3.45) 
Jo Jn Jo JdS 
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Combining (3.43), (3.44) with (3.45), we obtain the Eulerian derivative 



dJ2{n;V) = r f a 

Jo Jas 



^|curlyp + (D(y — g) ■ n) ■ (Dv ■ n — curly A n) 



Vr, dsdt. 



Finally we arrive at the expression (3.41) for the shape gradient VJ2. □ 

4 Gradient algorithm and numerical simulation 

In this section, we will give a gradient type algorithm and some numerical examples in 
two dimensions to prove that our previous methods could be very useful and efficient 
for the numerical implementation of the shape optimization problems for the unsteady 
Navier-Stokes flow. For the sake of simplicity, we only consider the cost functional 

4.1 A gradient type algorithm 

As we have just seen, the general form of the Eulerian derivative is 



dJ(0;y) = / / VJ- Fdsdt, 
Jo J as 



where V J denotes the shape gradient of the cost functional J. Ignoring regularization, 
a descent direction is found by defining 

V = -hkVJ (4.1) 

and then we can update the shape O as 

^k = {l + huV)^ (4.2) 

where hk is a descent step at fe-th iteration. 

There are also other choices for the definition of the descent direction. Since the 
gradient of the functional has necessarily less regularity than the parameter, an iterative 
scheme like the method of descent deteriortates the regularity of the optimized param- 
eter. We need to project or smooth the variation into H^{^1)'^. Hence, the method 
used in this paper is to change the scalar product with respect to which we compute a 
descent direction, for instance, ff^(J7)^. In this case, the descent direction is the unique 
element d £ if^(r2)^ such that at a fixed time t G [0,r] and for every V € H^{Q)'^, 



Bd:BVdx= VJ-V ds. (4.3) 

n Jas 

The computation of d can also be interpreted as a regularization of the shape gra- 
dient, and the choice of H^{^1)'^ as space of variations is more dictated by technical 
considerations rather than theoretical ones. 

The resulting algorithm can be summarized as follows: 
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(1) Choose an initial shape Oqi i-e-i choose an initial shape of dS since dD is fixed 
in our problem; 

(2) Compute the state system (2.4) and adjoint state system (3.36), then we can 
evaluate the descent direction by using (4.3) with = O^; 

(3) Set = (Id — /ifcdfc) J^fc, where /i^ is a small positive real number. 

The choice of the descent step hk is not an easy task. Too big, the algorithm is 
unstable; too small, the rate of convergence is insignificant. In order to refresh hk, we 
compare hk with hk^i- If {dk^dk-iju^ is negative, we should reduce the step; on the 
other hand, if dk and d^-i are very close, we increase the step. In addition, if reversed 
triangles are appeared when moving the mesh, we also need to reduce the step. 

In our algorithm, we do not choose any stopping criterion. A classical stopping 
criterion is to find that whether the shape gradients V J in some suitable norm is small 
enough. However, since we use the continuous shape gradients, it's hopeless for us to 
expect very small gradient norm because of numerical discretization errors. Instead, 
we fix the number of iterations. If it is too small, we can restart it with the previous 
final shape as the initial shape. 

4.2 Numerical examples 



To illustrate the theory, we want to solve the following minimization problem 

1 

/o Jn 



mjn^ / / \y-yd\^dxdt (4.4) 



subject to 

' dty - aAy + T>y ■y + Vp = f in x (0, 1) 
divy = infix (0,1) 

y = on 55 X (0, 1) (4.5) 

y = y^ ondD X (0, 1) 

j/(0) = inn. 
Where D := {{x,y) e : + < 0.64}, and the shape of the body 5 is to 
be optimized. We choose the velocity y^ = {O.lby, —0.15x)^ and the body force 
/ = (/i,/2)^: 

_ 453; atyilbx"^ + Iby"^ - 1) 

^ ~3UMT7 ^ 5(x2 + y2)3/2 



+ —t^x[ -46 - 25x2 - 252/2 _ + + 60 ^/x^ + ; 

25 \^ x^ + y^ + 2/2 ) 

45y atx(15x2 + 152/2 - 1) 



3iy^^T^ 5(x2 + 2/2)3/2 



lo/ 9 9 1 12 



+ -t2y -46 - 25x2 - 252/2 - ^— ^ + + 

25 \ x2 + 2/2 ^x2 + 2/2 
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The target velocity is determined by the data /, and the target shape of 
the domain Q. Our aim is to recover the shape of S which is a circle: dS = {{x,y) : 

a;2 + y2 = o.04}. 

The Navier-Stokes system (2.4) and the adjoint system (3.36) are discretized by 
using a mixed finite element method. Time discretization is effected using the backward 
Euler method and we assume that the time interval [0, 1] is divided into equal intervals 
of duration At = 0.05. Spatial discretization is effected using the Taylor-Hood pair 
[13] of finite element spaces on a triangular mesh, i.e., the finite element spaces are 
chosen to be continuous piecewise quadratic polynomials for the velocity and continuous 
piecewise linear polynomials for the pressure. Our numerical solutions are obtained 
under FreeFem++ [12] and we run the program on a home PC. 

We choose the initial shape of S to be elliptic: {{x,y) : + = 1/25}, and 

the initial finite element mesh was shown in Figure 4.1. 



Figure 4.1: Initial mesh with 125 nodes. 

Figure 4.2 — Figure 4.4 give the comparison between the target shape with iterated 
shape for the viscosity coefficients a = 0.1,0.01 and 0.001, respectively. In case of 
a = 0.1,0.01, we have fine results in Figure 4.2 and Figure 4.3. Unfortunately, we can 
not get a nice reconstruction for a = 0.001 as in Figure 4.4. 

Figure 4.5 represents the fast convergence of the cost functional for the various 
viscosity coefficients a = 0.1,0.01 and 0.001. 

5 Conclusion 

In this paper, the shape optimization in the two dimensional time-dependent Navier- 
Stokes flow has been presented. We employed the weak implicit function theorem to 
obtain the existence of the weak Piola material derivative, then we gave the descrip- 
tion of the shape derivative. Hence we derived the structures of shape gradients with 
respect to the shape of the variable domain for some time-dependent cost functionals 
by introducing the associated adjoint state system. A gradient type algorithm is effec- 
tively used for the minimization problem in various Reynolds number flows. Further 
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Figure 4.2: a = 0.1, CPU time: 124.531 s. 
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Figure 4.3: a = 0.01, CPU time: 120.125 s. 
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Figure 4.4: a = 0.001, CPU time: 622.813 s. 
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Figure 4.5: Convergence history for a = 0.1,0.01 and 0.001. 
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research is necessary on efficient implementations for very large Reynolds numbers and 
real problems in the industry. 
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